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ABSTRACT 

We use direct N-body simulations to study the inspiral and internal evolution 
of dense star clusters near the Galactic center. These clusters sink toward the 
center due to dynamical friction with the stellar background, and may go into 
core collapse before being disrupted by the Galactic tidal field. If a cluster 
reaches core collapse before disruption, its dense core, which has become rich in 
massive stars, survives to reach close to the Galactic center. When it eventually 
dissolves, the cluster deposits a disproportionate number of massive stars in the 
innermost parsec of the Galactic nucleus. Comparing the spatial distribution and 
kinematics of the massive stars with observations of IRS 16, a group of young He 
I stars near the Galactic center, we argue that this association may have formed 
in this way. 

Subject headings: Galaxy: center-Galaxy: nucleus-black hole physics — globular 
clusters: individual (Arches, Quintuplet) — stellar dynamics-methods: N-body 

1. Introduction 

Krabbe et al. (1995) found ~ 15 bright He I emission line stars within about 1 pc of the 
Galactic center, accompanied by many less luminous stars of spectral types O and B (Genzel 
et al. 2000). Genzel et al. (2000) have measured accurate positions and velocities of 41 early 
type stars in this region, and report proper motions for 26 of them. These stars are part 
of the co-moving group IRS 16, which was apparently formed 7-8 Myr ago in a starburst of 
mass ^ 10 4 M Q (Tamblyn & Rieke 1993). They show a high degree of anisotropy; most of 
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the He I stars in the Galactic center are on tangential orbits (Genzel et al. 2000). Detailed 
spectroscopic analysis of these Galactic center objects (Najarro et al. 1994) indicates that 
they are highly evolved, with a high surface ratio of helium to hydrogen uhc/^h — 1-1.67. 
Allen et al. (1990) classify them as Ofpe/WN9 stars, while Najarro et al. (1997) identify them 
as luminous blue variables (LBVs) with masses between 60 and 100 M Q . 

LBVs are the late evolutionary stages of very massive ( ^ 40 M Q ) stars (Langer et 
al. 1994). Massive stars remain in this stage for only a short while (~ 3 x 10 4 years) after 
leaving the main sequence and the helium-rich WN stage, placing these stars in a very narrow 
age bracket: 3.2-3.6 Myr (Langer et al. 1994). If these objects are lower-mass (25-40 M ) 
Wolf-Rayet stars, they may be somewhat older (5 to 7 Myr, see Testor et al. 1993), which is 
consistent with the age estimate of 7-8 Myr obtained from the model calculations of Tamblyn 
& Rieke (1993) and the independent age determination of 3-7 Myr by Krabbe et al. (1995). 
In either case, a firm age limit of ~ 7 Myr is indicated for IRS 16. The absence of detectable 
X-ray emission from these stars (Baganoff et al. 2001a, 2001b) argues in favor of the LBV 
interpretation, in which case the age limit drops to ~ 3.5 Myr. 

While the age of the IRS 16 group is fairly well constrained, the location at which 
it formed is not. One obvious possibility is that the starburst occurred at roughly the 
Galactocentric radius where the group is now observed. However, this model is problematic, 
as the formation of stars within a parsec of the Galactic Center is difficult. The tidal field of 
the central black hole is sufficient to unbind gas clouds with densities ^ 10 10 cm~ 3 (Morris 
1993). At a distance ^ 1 pc star formation is still easily prevented, even though the potential 
of the bulge starts to dominate over that of the black hole. 

Gerhard (2001) proposed that a massive star cluster of mass m formed at a distance 
of 30(m/10 6 M Q ) pc from the Galactic center can spiral in to the Galactic center by 
dynamical friction before being disrupted by the tidal field of the Galaxy or its own internal 
evolution. In order to survive in the Galactic central region the cluster core density has to 
exceed p c ^ 10 7 M Q pc~ 3 . It is unlikely that a star cluster would be born with such a high 
central density, but it may evolve into this state when core collapse occurs. However, even 
then the cluster must have been initially quite compact. Core collapse of a cluster boosts the 
central density, but can be strongly affected by mass loss from the cluster tidal boundary. 
In the strong tidal environment of the Galactic center, mass loss from the cluster perimeter 
may prevent core collapse altogether. 

We simulate dense star clusters using a direct TV-body approach, taking the external po- 
tential of the Galaxy and the effect of dynamical friction into account. Within this model we 
study the possibility that a cluster may go into core collapse before dynamical friction causes 
it to spiral in to the Galactic center. We include the dynamical friction term analytically, 
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applying it to the bound cluster mass (see Binney & Tremaine 1987). In §2 we discuss cluster 
evolution and dynamical friction, in order to interpret the results of our model calculations, 
which are presented in §3. We summarize in §4. 



2. Cluster dynamics 

2.1. Time scale for core collapse 

The dynamical evolution of a star cluster drives it toward core collapse (Antonov 1962; 
Spitzer & Hart, 1971) in which the central density runs away to a formally infinite value in 
a finite time. Core collapse occurs at 

^cc — C ^rlx j ( 1 ) 

where t r \ x is the cluster's "half-mass" relaxation time, 

rlx (GW) 1/2 8mA' [) 

Here G is the gravitational constant, n is the number of stars in the cluster, m is the total 
cluster mass, and r vir is the cluster's virial radius. The Coulomb logarithm In A ~ ln(O.ln) ~ 
10 typically. 

In an isolated cluster in which all stars have the same mass, c ~ 15 (Cohn 1980). In 
a multi-mass system, core collapse is determined by the accumulation of the most massive 
stars in the cluster center (Vishniac 1978; see also Chernoff & Weinberg 1990). We have 
performed direct iV-body simulations to determine the moment at which core collapse occurs, 
and hence the value of c in Eq. 1. 

The initial conditions of our model cluster are presented in Table 1. The cluster consists 
of 65536 stars distributed initially in a King (1966) model with King parameter Wo = 3. 
Each of the stars is randomly assigned a mass drawn from a Scalo (1986) initial mass function 
between 0.3 and 100 M , irrespective of position. The entire cluster is then rescaled to virial 
equilibrium. We choose a virial radius r vir = 0.167pc. These choices mimic the young 
Arches and Quintuplet star clusters, which are located somewhat farther (~ 30 pc) from the 
Galactic center. The resulting initial parameters (total mass, core radius, half mass radius, 
crossing time and relaxation time) are also listed in the table. 

Visual inspection of the core radius as a function of time indicates that core collapse 
occurs around t = 0.76 Myr, near the moment when the hard binary containing the most 
massive star reaches a binding energy of E < — 100 kT (where kT is the thermodynamic 
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Table 1: Initial conditions for the simulated star clusters. The columns give the number of 
stars, mass (in M ), the King parameter Wo, the core-, half mass-, virial- and tidal radii (all 
in parsecs). The last two columns give the initial half mass crossing time and the two-body 
relaxation (both in millions of years). 

n m W r corc r hm r vir r ti de t hm t r[x 

[M©] [pc] - [Myr] - 

65536 62549 3 0.116 0.140 0.167 0.525 0.0115 3.87 



energy scale of the stellar system; the total kinetic energy of the cluster is |n/cT). This 
binary was formed somewhat earlier (at t = 0.58 Myr), but at that time we could not 
identify the core as collapsed, as the core radius continued to contract. A little later (t = 
0.84 Myr), this binary is strongly perturbed by another star, resulting in a collision. Based 
on this information, we conclude that this particular simulation experienced core collapse 
at t ~ 0.76 Myr, so c ~ 0.20, which is consistent with Portegies Zwart & McMillan (2002), 
but somewhat larger than the c ~ 0.15 found in the Fokker-Planck simulations of Chernoff 
& Weinberg (1990). 



2.2. Inspiral to the Galactic center 

The mass M of the Galaxy within the cluster's orbit at distance R ( ^ 500 pc) from the 
Galactic center is taken to be (Sanders & Lowinger 1972; Mezger et al. 1996) 

M{R) = AR a , (3) 

where A = 4.25 x 10 6 M©(lpc)" a and a = 1.2. This slope fits the observed light distribution 
with constant M/L and the rotation curve derived from OH/IR stars and 21cm line obser- 
vations (Mezger et al. 1996). Earlier observations, however, claim a slightly shallower slope 
(Catchpole, Whitelock & Glass 1990). For clarity we adopt a = 1.2 for the remainder of this 
paper. The density at distance R then is 

. 1 dM Aa , 

Following Binney & Tremaine (1987), we find that the inspiral of the cluster towards the 
center due to dynamical friction is described by (McMillan & Portegies Zwart, 2002, ApJ in 
press) 

d -^ = -2\nA^-(^] 1/2 mR^\ (5) 



dt a + 1 \ A 
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Here m is constant and In A ~ hv(R/{r)) ~ 5 is a Coulomb logarithm (where (r) is the 
object's characteristic radius — roughly the half- mass radius in the case of a cluster) and 
X ~ 0.3 is a parameter which depends on the velocity of the cluster and the velocity dispersion 
of the stellar surroundings. In this case, %ln A ~ 1. The adopted value of In A is consistent 
with results from N-body simulations (Spinnato et al. 2003), who derive In A = 6.6±0.6 for a 
massive compact object that spirals-in, and somewhat smaller than the value In A ~ ln0.4iV 
used by Gerhard (2001), where N the number of stars with which the cluster interacts. For 
simplicity we write Eq. 5 as 

Rl^dR = - 7 di. (6) 



Solving the differential equation Eq. 5 with R(0) = Ri at time t — results in 

2/(3+a) 



R(t) = Ri 



(3 + aOT 



2R 



i(3+a) 



(7) 



Inverting this equation with R = Rf (the disruption radius) at t — t^f and substituting Eq. 3 
gives 

a + 1 1 _! /M(i2,-^ 1/2 



tdf — 



-m 



R! /2 - kR 



r] . 



a(a + 3) x In A 

where « = (R f /Ri) a/2 . For a = 1.2, xlnA = 1, and Rf <^ Ri, this becomes 



(8) 



. „ - ( m 

t df = 1.34 — 

V 10 4 M G 







- : / R . \ (3+a)/2 



lpc 



Myr . 



(9) 



3. Results 

In order to test the hypothesis that a cluster can experience core collapse before reaching 
the Galactic center, it is instructive to compare the dynamical friction inspiral time scale 
with the time scale for internal cluster evolution. We define 



t, 



cc 

a{a + 3) cxlnA / m\ V2 m r vir 3 / 2 

" a + 1 81nA \MJ (m) R 3 / 2 _ K tfl 2 ' 1 ' 

1 j 

Here (m) is the mean mass of cluster stars. For small Rf, this reduces to 

/0.29c X hiA\ /my/2 /r vir \ 3/2 
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There are now three distinct regimes: (i) If 77 •< 1 (far from the Galactic center: 
[i?/pc] 4 - 2 > [n/2.1 x 10 3 ] 2 [m/M ] [r vir /pc] 3 ), the cluster core collapses essentially at its 
original distance from the Galactic center; thereafter it dissolves, mainly by tidal stripping 
and mass loss due to stellar evolution, at constant Galactocentric radius, (ii) If 77 ~ 1 (inter- 
mediate distance to the Galactic center), cluster inspiral and core collapse occur on about 
the same time scale, (iii) If 77 ^> 1 (close to the Galactic center), the cluster spirals in without 
significant internal evolution. 

For example, substituting the initial conditions of Tab. 1 (m = 65 000M Q , r vir = 
0.167pc) and Eq. 3 into Eq. 11, we can write 

77 ~ 150c X !^iT 21 . (12) 
In A 

Taking c^lnA/lnA ~ 0.13, we find that this cluster will experience core collapse before it 
reaches the Galactic center if it was born at Ri 4pc. 

More generally, Fig. 1 shows, as functions of Galactocentric radius, the virial radii (solid 
curve) and an estimate for the initial tidal radii (dots) of star clusters with masses and 
structure parameters as listed in Tab. 1, for clusters with 77 = 1 in Eq. 11 (core collapse upon 
arrival at the Galactic center), and take c^lnA/lnA ~ 0.13. The dashed curve presents 
an estimate of the Jacobi radius of the star cluster in the Galactic tidal field (see Eq. 4 in 
Gebhard 2001, or Eq. 24 in McMillan & Portegies Zwart 2002). A 65 000M© star cluster 
which is born with parameters to the right or below the solid curve is expected to experience 
core collapse before it reaches the Galactic center. 

The circles and bullets at r vir = 0.167pc in Fig. 1 indicate the outcomes of simulations 
performed with cluster initial conditions as presented in Tab. 1, but with varying values of 
initial Galactocentric radius Ri. The initial models were placed at Ri = 2, 3, 4, 5, 6, and 
8pc. Bullets indicate that a model experienced core collapse before reaching the Galactic 
center; circles indicate disruption before significant contraction of the cluster core. 

The equations of motion of the 65536 (64k) stars in the simulations were computed us- 
ing the Starlab software environment which combines the iV-body integrator kira and the 
binary evolution package SeBa (Portegies Zwart et al. 2001) 5 . The Galactic tidal field and 
the effects of dynamical friction were taken properly into account by solving Eq. 5 numeri- 
cally during the integration of the equations of motion. In these calculations, the clusters 
could lose mass by tidal stripping, high velocity stellar ejections, and stellar winds. At any 
moment in time we determined the total cluster mass from all bound stars; this may slightly 



5 see http://manybody.org 
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Fig. I. — Virial radius r vir of a cluster for which core collapse is expected by the time of 
arrival in the Galactic center (solid curve). Here the mass of the cluster is 65OOOM , with 
64k stars. The dotted curve gives 3r vir which corresponds to the radius where the cluster 
density goes to zero for a King model with Wo = 3. The dashed line gives the Jacobi 
radius of the cluster in the tidal field of the Galaxy. Circles (lower left) and bullets show the 
simulations performed for this study. Circles indicate that the cluster dissolved before core 
collapse; bullets indicate simulations which did experience core collapse. 
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underestimate the friction force. The dynamical friction term was then applied to each of 
the bound stars, but not to unbound stars. The calculations were carried out using the 
special-purpose GRAPE-6 computer (Makino et al. 1997; Makino 2001). 



Table 2: 



Ri 


^cc 


Rcc 


^diss 


Rf 


[PC] 


[Myr] 


M 


[Myr] 


[PC] 


2 






1.08 


1.1 


3 






1.01 


1.1 


4 






0.83 


1.3 


5 


0.65 


3.3 


1.19 


1.8 


6 


0.68 


4.6 


> 1.44 


< 2.5 


8 


0.56 


7.0 


> 1.05 


< 6.4 



The results of the simulations in Tab. 2, may be summarized as follows. The models 
with Ri = 2 and 3pc both dissolved at Rf = 1.1 pc. They did not experience core collapse, 
nor were any persistent hard binaries formed. (For defmiteness, we take a cluster to have 
dissolved once the bound mass drops below 6 000M Q .) The model with Ri = 4pc dissolved 
at Rf = 1.3 pc, but core collapse in this case is uncertain. A few hard (E < — 10 kT) 
binaries formed at t = 0.51 Myr, at a distance of R = 2.3 pc. One of these binaries hardened 
to E ^ — 50 kT at t = 0.63 Myr and R = 1.8 pc; the cluster dissolved a little later, at 
t = 0.83 Myr. This ambiguity is consistent with the cluster's location close to the solid curve 
in Fig. 1. The models with Ri > 5pc all experienced core collapse at t cc ~ 0.6 Myr. The 
Ri = 5pc model dissolved at t = 1.19 Myr at a galactocentric distance of Rf = 1.8 pc. 6 . The 
other models were not continued to the point of dissolution. 

Figure 2 shows a top view of the orbit of the cluster with Ri = 5 pc. The dotted curve 
indicates the expected orbit of the cluster if its mass would remain constant and was not 
affected by stellar evolution, internal relaxation or by the external tidal field. This constant- 
mass point spirals-in slightly more quickly than the cluster in which mass loss is taken into 
account self consistently (solid curve). 

Figure 3 shows four subsequent snapshots (gray shadings and contours) at time intervals 
of 0.4 Myr for the cluster with R { = 5 pc. 



6 An animation of this simulation is available at http://manybody.org/starlab.html. See also 
http : / /www. ids . ias . edu/^starlab/animations/. 
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Fig. 2. — Orbital evolution of the star cluster with Ri = 5pc (solid curve). The cluster is 
tracked via its density center, which becomes hard to determine accurately when the cluster 
contains only a few stars, i.e: near the Galactic center. The star indicated the moment of 
core collapse. The dotted curve shows the orbit for a constant-mass point with the same 
mass as the cluster at birth. The bullets in the orbit indicate 10 5 year time intervals for the 
cluster orbit. In the dotted curve these intervals are indicated with small circles. 
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Fig. 3. — Top view of model Ri = 5pc at t = 0.1 Myr, 0.4 Myr, 0.8 Myr and 1.2 Myr. Density 
is gray shaded linearly between maximum (dark) and zero density (light) scaled individually 
to each panel. The contours indicate a constant stellar density of 10 stars/pc 3 , 50 stars/pc 3 , 
100 stars/pc 3 , etc 7 . (Note that this calculation was performed with a different value of In A 
than was adopted in Fig. 2.) 
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lt is at first somewhat surprising that Eq. 11 agrees so well with the simulation results, 
as the cluster mass in the latter is a function of time which is neglected in the analytic 
form. Gerhard (2001) corrects for mass loss using an isothermal model which implies that 
the cluster loses mass at a constant rate until disruption. This would result in a factor 
of two increase in the dynamical friction time. It turns out that this overestimates mass 
loss considerably (see fig. 2). Most mass in the Wo=3 King model is lost near the end of 
the cluster lifetime when the stellar density in the environment becomes comparable to the 
cluster density near the half mass radius. 

In the clusters which did not experience core collapse, stars at disruption were spread 
over a broad range in radii: Rf £ R £ Ri pc. Stars more massive than 40 M Q were not 
distributed in a significantly different way from low-mass stars. However, in the Ri = 5pc 
model, which did reach core collapse before disruption, the massive stars became much more 
centrally concentrated than the other cluster members. The more massive stars penetrated 
closer to the Galactic center because they sank to the cluster core, whereas low-mass stars 
were lost from the cluster at an earlier stage, when the cluster was farther from the center. 

While we have shown that core collapsed clusters preferentially deposit their most mas- 
sive stars closest to the Galactic center, our models differ in some important ways from the 
observations reported by Genzel et al. (2000). The end-point of our Ri = 5 core-collapse 
cluster is a tight clump of stars with spatial extent (half-mass radius) comparable to those 
observed, but the clump was deposited at a radius of almost 2 pc, not within the inner- 
most 1 pc, as is usually assumed (Krabbe et al. 1995). However, we expect that more 
massive clusters (starting with smaller virial radii or larger galactocentric radii — see Eq. 11), 
more centrally concentrated systems or systems on elliptic orbits, penetrate deeper into the 
Galactic potential. 

Genzel et al. (2000) present (their Fig. 5) the velocity anisotropy of 12 of the He I stars 
near the Galactic center. The anisotropy parameter is 7 = (i>f — v^)/(vf + v%), where v t 
and v r are the transverse and radial velocity components of the stellar proper motions. The 
average velocity anisotropy of these stars is (7) = 0.59 ± 0.48. When we exclude the one 
star with an unusually low value of 7 = —0.83 ± 0.33, the average anisotropy increases to 
(7) = 0.72 ±0.18. We measured the velocity anisotropy among the star with a mass > 40 M Q 
of our Ri = 5pc model at the moment it disrupted, and found (7) = 0.79 ± 0.23. For all 
stars in the cluster the mean velocity anisotropy is (7) = 0.04 ± 0.63, which is consistent 
with being isotropic. Likewise the sky-projected radial and tangential velocities of all 104 
proper motion stars in the sample of Genzel et al. (2000) is consistent with overall isotropy. 

If IRS 16 is indeed a remnant cluster core, our simulations provide no easy explanation 
of the rather broad stellar distribution perpendicular to the supposed cluster orbit plane, 
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nor for the large dispersion in their velocities. The stars in our simulations are eventually 
spread out in the orbit plane, they are quite tightly confined in the direction perpendicular 
to this plane. The dispersion in the velocity distribution of these stars then would be of 
the order of the cluster velocity dispersion: of the order of ten kilometers per second, rather 
than the observed dispersion of a few lOOkm/s. The disruption of two star clusters in 
short succession would not reproduce all kinematic information. We speculate that other 
dynamical processes, such as the effects of primordial binaries, the presence of a central black- 
hole binary or possible inhomogeneities in the background potential, might have operated in 
IRS 16 to increase its scale height out of the plane and to drive the velocity dispersion to its 
observed values. At present, however, we have no ready solution to this conundrum. 



4. Summary 

We have critically examined the hypothesis proposed by Gerhard (2001) that the group 
IRS 16 may be the remnant of a much larger cluster that formed farther from the Galactic 
center and sank toward the center via dynamical friction. IRS 16 contains about 40 early- 
type stars, including at least 15 very luminous He I stars, all lying within ~ 0.4 pc of the 
Galactic center. 

We have studied this possibility by performing a series of direct iV-body simulations in 
which dynamical friction is taken into account in a semi-analytic fashion, and included self- 
consistently in the equations of motion of the cluster stars. The iV-body calculations were 
performed with 65536 stars and were run on the GRAPE-6. Stellar masses were selected 
from a realistic mass function and stars were initially distributed as a Wo = 3 King model 
with virial radius r vir = 0.167pc. We find that, in order for a clump of massive stars to 
survive, the cluster must have experienced core collapse during the inspiral. Core collapse 
deposits the observed high proportion of early type stars close to the Galactic center, and 
prevents a spread of massive stars to larger distances. The anisotropy observed for the early- 
type stars in IRS 16 is consistent with our model calculations. However, the spatial extent 
and high dispersion in the velocities of the observed cluster are not satisfactorily explained 
with the current simulations. The presence of primordial binaries or a binary of intermediate 
mass black holes in the cluster center may be required to explain these observations. 

Our approximation to the dynamical friction term has some limitations, as the parame- 
ter xhiA is fixed in our simulations. In reality, this term may differ from the adopted value, 
will probably vary with time, and may depend on a number of external factors. More accu- 
rate measurements of this parameter are presented by Spinnato et al. (2003). Regardless 
of this uncertainty, we are still able to draw some firm conclusions. We find that dense star 
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clusters in a strong tidal field experience core collapse on a time scale similar to that for an 
isolated cluster, and that core collapse can occur before a cluster near the Galactic center is 
tidally disrupted. 

When the cluster experiences core collapse, the fraction of massive stars deposited near 
the center is much greater than when core collapse is averted by tidal disruption. Our 
simulated clusters dissolve when their core densities fall below a few million M pc~ 3 , more 
than an order of magnitude higher than the local background stellar density. Our calculations 
were performed using rather low concentration Wo=3 clusters. We expect that clusters with 
higher initial concentrations would penetrate deeper and more easily to the Galactic center. 
Variations in the initial orbit of the cluster may also prove to be efficient in transporting 
stars closer to the Galactic center. Note that our choice of initial conditions are possibly 
among the least favorable to explain the observations. The parameter space for clusters 
which experience core collapse before reaching the Galactic center may therefore be even 
larger than suggested here. 

Our model calculations support the scenario proposed by Gerhard (2001) to explain the 
presence of a population of early type stars within a parsec of the Galactic center. If born 
at a distance of ~ 5pc, the primordial cloud from which the cluster formed should have had 
an initial density on the order of 10 6 cm -3 , but this density might be lower if the stars in 
IRS 16 originated in a somewhat more massive cluster at a greater distance. 
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